!  http://www.cs.rpi.edu/~szymansk/OOF90/F90_Objects.html
module lubricant_class
    implicit none
    integer, parameter :: viscous = 1
    integer, parameter :: maxwellviscoelastic = 2
    integer, parameter :: newtonian = 1
    integer, parameter :: rheeeyering = 2
    integer, parameter :: bairviscoplastic = 3
    integer, parameter :: barusdirecttemperature = 1
    integer, parameter :: barusinversetemperature = 2
    integer, parameter :: roelands = 3
    integer, parameter :: linearpressuredirecttemperature = 1
    integer, parameter :: linearpressureinversetemperature = 2
    integer, parameter :: exponentialpressuredirecttemperature = 3
    integer, parameter :: exponentialpressureinversetemperature = 4
    logical, parameter :: OK = .true.
    logical, parameter :: notOK = .false.

    type lub_t
        integer :: rheologybehaviour
        integer :: highshearviscositybehaviour
        integer :: lowshearviscositybehaviour
        integer :: limitingshearstressdependence
        integer :: elasticshearmodulusdependence
        double precision :: lowshearviscosityT0, lowshearviscosity0, lowshearviscosityalphaZ, lowshearviscositybetaS
        double precision :: limitingshearstressT0, limitingshearstress0, limitingshearstressalpha, limitingshearstressbeta
        double precision :: elasticshearmodulusT0, elasticshearmodulus0, elasticshearmodulusalpha, elasticshearmodulusbeta
    end type lub_t




    contains




        elemental function highshearviscosity(lub, T, p, fxz, fxy, fyz, fxx, fyy, fzz)
            implicit none
            double precision :: highshearviscosity
            type (lub_t), intent (in) :: lub
            double precision, intent (in) :: T, p
            double precision, intent(in) :: fxz
            double precision, intent(in), optional :: fxy, fyz, fxx, fyy, fzz
            double precision :: gammaR, tauE, eta
            if (lub%highshearviscositybehaviour == newtonian) then
                highshearviscosity = lowshearviscosity(lub, T, p)
            elseif (lub%highshearviscositybehaviour == rheeeyering) then
                gammaR = seconddeviatoricinvariant(fxz, fxy, fyz, fxx, fyy, fzz)
                tauE = limitingshearstress(lub, T, p)
                eta = lowshearviscosity(lub, T, p)
                highshearviscosity = tauE / gammaR * asinh( eta * gammaR / tauE)
            elseif (lub%highshearviscositybehaviour == bairviscoplastic) then
                gammaR = seconddeviatoricinvariant(fxz, fxy, fyz, fxx, fyy, fzz)
                tauE = limitingshearstress(lub, T, p)
                eta = lowshearviscosity(lub, T, p)
                highshearviscosity = tauE / gammaR * (1 - exp(-gammaR*eta/tauE))
            else
                highshearviscosity = -1.d0
            end if
            return
        end function highshearviscosity




!        elemental function highshearviscosity(lub, T, p, fxz, fxy, fyz, fxx, fyy, fzz)
!            implicit none
!            double precision :: highshearviscosity
!            type (lub_t), intent (in) :: lub
!            double precision, intent (in) :: T, p
!            double precision, intent(in) :: fxz
!            double precision, intent(in), optional :: fxy, fyz, fxx, fyy, fzz
!            double precision :: tauR, tauE, eta
!            if (lub%highshearviscositybehaviour == newtonian) then
!                highshearviscosity = lowshearviscosity(lub, T, p)
!            elseif (lub%highshearviscositybehaviour == rheeeyering) then
!                tauR = seconddeviatoricinvariant(fxz, fxy, fyz, fxx, fyy, fzz)
!                tauE = limitingshearstress(lub, T, p)
!                eta = lowshearviscosity(lub, T, p)
!                highshearviscosity = eta * tauR / tauE / sinh(tauR/tauE)
!            elseif (lub%highshearviscositybehaviour == bairviscoplastic) then
!                tauR = seconddeviatoricinvariant(fxz, fxy, fyz, fxx, fyy, fzz)
!                tauE = limitingshearstress(lub, T, p)
!                eta = lowshearviscosity(lub, T, p)
!                highshearviscosity = -eta * tauR / tauE / log(1-tauR/tauE)
!            else
!                highshearviscosity = -1.d0
!            end if
!            return
!        end function highshearviscosity




        elemental function lowshearviscosity(lub, T, p)
            implicit none
            double precision :: lowshearviscosity
            type (lub_t), intent (in) :: lub
            double precision, intent (in) :: T, p
            double precision :: t0, x0, a, b
            if (lub%lowshearviscositybehaviour == barusdirecttemperature) then
                t0 = lub%lowshearviscosityT0
                x0 = lub%lowshearviscosity0
                a = lub%lowshearviscosityalphaZ
                b = lub%lowshearviscositybetaS
                lowshearviscosity = x0 * exp(a * p - b * (T - t0))
            elseif (lub%lowshearviscositybehaviour == barusinversetemperature) then
                t0 = lub%lowshearviscosityT0
                x0 = lub%lowshearviscosity0
                a = lub%lowshearviscosityalphaZ
                b = lub%lowshearviscositybetaS
                lowshearviscosity = x0 * exp(a * p + b * (1/T - 1/t0))
            elseif (lub%lowshearviscositybehaviour == roelands) then
                t0 = lub%lowshearviscosityT0
                x0 = lub%lowshearviscosity0
                a = lub%lowshearviscosityalphaZ
                b = lub%lowshearviscositybetaS
                lowshearviscosity = 6.31d-5 * exp( log(x0/6.31d-5) * (p/196.d6 + 1)**a * ((t0-138.15)/(T-138.15))**b )
            else
                lowshearviscosity = -1.d0
            end if
            return
        end function lowshearviscosity




        elemental function tangentpressureviscositycoefficient(lub, T, p)
            implicit none
            double precision :: tangentpressureviscositycoefficient
            type (lub_t), intent (in) :: lub
            double precision, intent (in) :: T
            double precision, intent (in), optional :: p
            double precision :: pp
            pp = 0.196d9
            if (present(p)) pp = p
            tangentpressureviscositycoefficient = log(lowshearviscosity(lub, T, pp) / lowshearviscosity(lub, T, 0.d0))/pp
            return
        end function tangentpressureviscositycoefficient




        elemental function limitingshearstress(lub, T, p)
            implicit none
            double precision :: limitingshearstress
            type (lub_t), intent (in) :: lub
            double precision, intent (in) :: T, p
            double precision :: t0, x0, a, b
            integer :: dependence
            dependence = lub%limitingshearstressdependence
            t0 = lub%limitingshearstressT0
            x0 = lub%limitingshearstress0
            a = lub%limitingshearstressalpha
            b = lub%limitingshearstressbeta
            limitingshearstress = anonymousparam(T, p, dependence, x0, t0, a, b)
            return
        end function limitingshearstress




        elemental function elasticshearmodulus(lub, T, p)
            implicit none
            double precision :: elasticshearmodulus
            type (lub_t), intent (in) :: lub
            double precision, intent (in) :: T, p
            integer :: dependence
            double precision :: t0, x0, a, b
            dependence = lub%elasticshearmodulusdependence
            t0 = lub%elasticshearmodulusT0
            x0 = lub%elasticshearmodulus0
            a = lub%elasticshearmodulusalpha
            b = lub%elasticshearmodulusbeta
            elasticshearmodulus = anonymousparam(T, p, dependence, x0, t0, a, b)
            return
        end function elasticshearmodulus




        elemental function anonymousparam(T, p, dependence, x0, t0, a, b)
            implicit none
            double precision :: anonymousparam
            double precision, intent (in) :: T, p
            integer, intent(in) :: dependence
            double precision, intent(in) :: t0, x0, a, b
            double precision :: x
            if (dependence == linearpressuredirecttemperature) then
                x = (x0 + a*p) * exp(-b*(T-t0))
                anonymousparam = max(0.d0, x)
            elseif (dependence == linearpressureinversetemperature) then
                x = (x0 + a*p) * exp(b*(1/T-1/t0))
                anonymousparam = max(0.d0, x)
            elseif (dependence == exponentialpressuredirecttemperature) then
                anonymousparam = x0 * exp(a*p - b*(T-t0))
            elseif (dependence == exponentialpressureinversetemperature) then
                anonymousparam = x0 * exp(a*p + b*(1/T-1/t0))
            else
                anonymousparam = -1.d0
            end if
            return
        end function anonymousparam




        elemental function seconddeviatoricinvariant(fxz, fxy, fyz, fxx, fyy, fzz)
            implicit none
            double precision :: seconddeviatoricinvariant
            double precision, intent(in) :: fxz
            double precision, intent(in), optional :: fxy, fyz, fxx, fyy, fzz
            double precision :: sxz, sxy, syz, sxx, syy, szz
            sxz = fxz
            sxy = 0.d0
            syz = 0.d0
            sxx = 0.d0
            syy = 0.d0
            szz = 0.d0
            if (present(fxy)) sxy = fxy
            if (present(fyz)) syz = fyz
            if (present(fxx)) sxx = fxx
            if (present(fyy)) syy = fyy
            if (present(fzz)) szz = fzz
            seconddeviatoricinvariant = sqrt( ((sxx-syy)**2.d0 + (syy-szz)**2.d0 + (szz-sxx)**2.d0)/6.d0 &
                + sxy*sxy + syz*syz + sxz*sxz )
            return
        end function seconddeviatoricinvariant




end module lubricant_class
